function [DGrid,VpGrid,VsGrid] = runsim_AVO_fixed_dens(simType, projectPath, inputPath,...
    outputPath, noIt,noSim,...
    vp_hd,vs_hd,fileid, bounds_den,bounds_vp,... % eliminated DenHD FT
    bounds_vs,XX,YY,ZZ,joint_VpVs,joint_VsD,krig_ANG_den,krig_RANGE_den,var_ANG_den,var_RANGE_den,...
    krig_ANG_vp,krig_RANGE_vp,var_ANG_vp,var_RANGE_vp,...
    krig_ANG_vs,krig_RANGE_vs,var_ANG_vs,var_RANGE_vs, varType, varNugget,...
    zoneFileName, nbrZones, nbrClassesBihist, covTab)
%% EXECUTE EXTERNALLY THE SIMULATIONS - DSS AND DSS-BISEC - RUN ONE SIM AT TIME
%
% noIt - NUMBER OF ITERATION
% noSim - NUMBER OF SIMULATION
% xMax - MAXIMUM GRID SIZE XX DIRECTION
% yMax - MAXIMUM GRID SIZE YY DIRECTION
% zMax - MAXIMUM GRID SIZE ZZ DIRECTION
% den_hd [string] - DENSITY HARD-DATA
% vp_hd [string] - VP HARD-DATA
% vs_hd [string] - VS HARD-DATA
% D_Vp_corr - GLOBAL CORRELATION BTW. DENSITY AND Vp
% Vp_Vs_corr - GLOBAL CORRELATION BTW. DENSITY AND VS
% fileid - LOG FILE ID
% bounds_den - MINIMUM AND MAXIMUM VALUES FOR DENSITY
% bounds_vp - MINIMUM AND MAXIMUM VALUES FOR Vp
% bounds_vs - MINIMUM AND MAXIMUM VALUES FOR Vs
% XX - GRID PARAMETERS FOR XX DIRECTION
% YY - GRID PARAMETERS FOR YY DIRECTION
% ZZ - GRID PARAEMTERS FOR ZZ DIRECTION
% joint_d_vp [string] - JOINT PROBABILITY DISTRIBUTION FILE DENSITY vs. Vp
% joint_vp_vs [string] - JOINT PROBABILITY DISTRIBUTION FILE Vp vs. Vs
%
% D_DSS [N x noSim] - ARRAY WITH ALL SIMULATIONS OF DENSITY
% Vp_coDSS [N x noSim] - ARRAY WITH ALL SIMULATIONS OF Vp
% Vs_coDSS [N x noSim] - ARRAY WITH ALL SIMULATIONS OF Vs
%
% FUNCTIONS:
% par_file_AVO.m - WRITE PARAMETER FILE FOR SIMULATION
%
% MODIFICATIONS:
% 04 APR 2014 - MAJOR PERFORMANCE UPDATE
% 25 NOV 2013 - VARIOGRAM WITH 2 STRUCT ARE NOW INPUT WITH CASE
% 16 JUN 2013 - NEW SIMULATION ALGORITHM - DSS PARALLEL C. NEW INPUTS FOR
% KRIGING AND VARIOGRAM PARAMETERS
% 12 MAR 2013 - NEW INPUTS FOR HARD-DATA FILES
% 31 OCT 2012 - NEW INPUTS FOR GRID PARAMETERS
% 30 OCT 2012 - NEW INPUTS FOR PROPERTY BOUNDS
% 04 SEP 2012 - CHANGES FOR SINGLE SIMULATIONS RUN - MEMORY UPDATE.
% 20 AUG 2012 - PERFORMANCE UPDATE
%
% LA - 21 AUG 2012
%

%% RUN INITAL SIMULATIONS

if noIt==1
    if (strcmp(simType ,'SGEMS') == 1)
        tic
        fprintf(fileid, '%s %d \n', 'Simulating VP - Simulation number: ', noSim); % WRITE NUMBER OF SIM
        VpGrid = gslib2grid(parFileDSS(simType, 1, 1, 0, [projectPath, inputPath], vp_hd, 1, ...
            [projectPath, outputPath, 'vp_DSS_',num2str(noIt),'_',num2str(noSim)], 0,...
            'No File', repmat('No File', nbrZones,1), bounds_vp,XX,YY,ZZ,'No File',0,'No File',krig_ANG_vp,krig_RANGE_vp,var_ANG_vp,var_RANGE_vp, varType, varNugget(:,1),...
            [projectPath, inputPath, zoneFileName], nbrZones,...
            nbrClassesBihist, 0, 'No File', 0, 'No file', 0, covTab),XX(1,1), YY(1,1), ZZ(1,1));

        fprintf(fileid, '%s %d \n', 'Simulating Vs - Simulation number: ', noSim); % WRITE NUMBER OF SIM
        VsGrid = gslib2grid(parFileDSS(simType, 1, 1, 1, [projectPath, inputPath], vs_hd, 1, [projectPath, outputPath, 'vs_coDSS_',num2str(noIt),'_',num2str(noSim)], 0,...
            'No File', joint_VpVs, bounds_vs,XX,YY,ZZ,'No File',0,[projectPath, outputPath,'vp_DSS_',num2str(noIt),'_',num2str(noSim),'_1.sgems'],...
            krig_ANG_vs,krig_RANGE_vs,var_ANG_vs,var_RANGE_vs, varType, varNugget(:,2), [projectPath, inputPath, zoneFileName], nbrZones,...
            nbrClassesBihist, 0, 'No File', 0, 'No file', 0, covTab),XX(1,1), YY(1,1), ZZ(1,1));     

        fprintf(fileid, '%s %d \n', 'USE LFM AS HD (FT)Simulating DEN - Simulation number: ', noSim); % WRITE NUMBER OF SIM
%        DGrid = gslib2grid(parFileDSS(simType, 1, 1, 1, [projectPath, inputPath], den_hd, 1, [projectPath,outputPath, 'd_coDSS_',num2str(noIt),'_',num2str(noSim)],0,...
%            'No File',joint_VsD, bounds_den,XX,YY,ZZ,'No File',0.9,[projectPath,outputPath,'vs_coDSS_',num2str(noIt),'_',num2str(noSim),'_1.sgems'],...
%            krig_ANG_den,krig_RANGE_den,var_ANG_den,var_RANGE_den, varType, varNugget(:,3), [projectPath, inputPath, zoneFileName], nbrZones,...
%            nbrClassesBihist, 1, 'No File', 0, 'No file', 0, covTab),XX(1,1), YY(1,1), ZZ(1,1));

        % FT load LFM density as simulation  matrix
        [~,dgrid] = import_gslib([projectPath 'input\DeFlat.gslib']);
        DGrid = gslib2grid(dgrid,XX(1,1), YY(1,1), ZZ(1,1));
        
            twtsim = toc;
        fprintf(fileid, '%s %d \n', 'Simulations - elapsed time (sec)', twtsim); % WRITE NUMBER OF SIM
        clear twtsim        
    else
%         tic
%         fprintf(fileid, '%s %d \n', 'Simulating Density - Simulation number: ', noSim); % WRITE NUMBER OF SIM
%         DGrid = gslib2grid(parFileDSS(simType, 1, 1, 0, [projectPath, inputPath], den_hd, 1, [projectPath,outputPath, 'd_DSS_',num2str(noIt),'_',num2str(noSim)], 0,...
%             'No File', repmat('No File', nbrZones,1), bounds_den,XX,YY,ZZ,'No File',0,'No File',krig_ANG_den,krig_RANGE_den,var_ANG_den,var_RANGE_den, varType,...
%             varNugget(:,1), [projectPath, inputPath, zoneFileName], nbrZones,...
%             nbrClassesBihist, 0, 'No File', 0, 'No file', 0, covTab),XX(1,1), YY(1,1), ZZ(1,1));
% 
%         fprintf(fileid, '%s %d \n', 'Simulating Vp - Simulation number: ', noSim); % WRITE NUMBER OF SIM
%         VpGrid = gslib2grid(parFileDSS(simType, 1, 1, 1, [projectPath, inputPath], vp_hd, 1, [projectPath,outputPath, 'vp_coDSS_',num2str(noIt),'_',num2str(noSim)], 0,...
%             'No File', joint_d_vp, bounds_vp,XX,YY,ZZ,'No File',0,[projectPath,outputPath,'d_DSS_',num2str(noIt),'_',num2str(noSim),'_1.out'],...
%             krig_ANG_vp,krig_RANGE_vp,var_ANG_vp,var_RANGE_vp, varType, varNugget(:,2), [projectPath, inputPath, zoneFileName], nbrZones,...
%             nbrClassesBihist, 0, 'No File', 0, 'No file', 0, covTab),XX(1,1), YY(1,1), ZZ(1,1));     
% 
%         fprintf(fileid, '%s %d \n', 'Simulating Vs - Simulation number: ', noSim); % WRITE NUMBER OF SIM
%         VsGrid = gslib2grid(parFileDSS(simType, 1, 1, 1, [projectPath, inputPath], vs_hd, 1, [projectPath,outputPath, 'vs_coDSS_',num2str(noIt),'_',num2str(noSim)],0,...
%             'No File', joint_VsD, bounds_vs,XX,YY,ZZ,'No File',0,[projectPath,outputPath,'vp_coDSS_',num2str(noIt),'_',num2str(noSim),'_1.out'],...
%             krig_ANG_vs,krig_RANGE_vs,var_ANG_vs,var_RANGE_vs, varType, varNugget(:,3), [projectPath, inputPath, zoneFileName], nbrZones,...
%             nbrClassesBihist, 0, 'No File', 0, 'No file', 0, covTab),XX(1,1), YY(1,1), ZZ(1,1));
% 
%         twtsim = toc;
%         fprintf(fileid, '%s %d \n', 'Simulations - elapsed time (sec)', twtsim); % WRITE NUMBER OF SIM
%         clear twtsim                
%         
    end
%%  RUN j=1++ SIMULATIONS    
else  
    if (strcmp(simType ,'SGEMS') == 1);
        tic
        fprintf(fileid, '%s %d \n', 'Simulating VP - Simulation number: ', noSim); % WRITE NUMBER OF SIM
        VpGrid = gslib2grid(parFileDSS(simType, 1, 1, 0, [projectPath, inputPath], vp_hd, 1, [projectPath,outputPath, 'vp_DSS_',num2str(noIt),'_',num2str(noSim)], 5,...
           [ projectPath,outputPath, 'best_vp_',num2str(noIt-1),'.sgems'], repmat('No File', nbrZones,1), bounds_vp,XX,YY,ZZ, ...
           [projectPath,outputPath, 'best_cc_vp',num2str(noIt-1),'.sgems'],0,'No File',krig_ANG_vp,krig_RANGE_vp,var_ANG_vp,var_RANGE_vp, varType, varNugget(:,1),...
           [projectPath, inputPath, zoneFileName], nbrZones,...
            nbrClassesBihist, 0, 'No File', 0, 'No file', 0, covTab),XX(1,1), YY(1,1), ZZ(1,1));

        fprintf(fileid, '%s %d \n', 'Simulating Vs - Simulation number: ', noSim); % WRITE NUMBER OF SIM       
        VsGrid = gslib2grid(parFileDSS(simType, 1, 1, 1, [projectPath, inputPath], vs_hd, 1, [projectPath,outputPath, 'vs_coDSS_',num2str(noIt),'_',num2str(noSim)], 5,...
           [ projectPath,outputPath, 'best_vs_',num2str(noIt-1),'.sgems'], joint_VpVs, bounds_vs,XX,YY,ZZ, ...
           [projectPath,outputPath, 'best_cc_vs',num2str(noIt-1),'.sgems'],0,[projectPath,outputPath, 'vp_DSS_',num2str(noIt),'_',num2str(noSim),'_1.sgems']...
           ,krig_ANG_vs,krig_RANGE_vs,var_ANG_vs,var_RANGE_vs, varType, varNugget(:,2), [projectPath, inputPath, zoneFileName], nbrZones,...
            nbrClassesBihist, 0, 'No File', 0, 'No file', 0, covTab),XX(1,1), YY(1,1), ZZ(1,1));

        fprintf(fileid, '%s %d \n', 'USE LFM AS HD (FT)Simulating DEN - Simulation number: ', noSim); % WRITE NUMBER OF SIM
%        DGrid = gslib2grid(parFileDSS(simType, 1, 1, 1, [projectPath, inputPath], den_hd, 1, [projectPath,outputPath, 'd_coDSS_',num2str(noIt),'_',num2str(noSim)], 5,...
%           [projectPath, outputPath, 'best_D_',num2str(noIt-1),'.sgems'], joint_VsD , bounds_den,XX,YY,ZZ, ...
%           [projectPath,outputPath, 'best_cc_den',num2str(noIt-1),'.sgems'],0.9,[projectPath,outputPath, 'vs_coDSS_',num2str(noIt),'_',num2str(noSim),'_1.sgems']...
%           ,krig_ANG_den,krig_RANGE_den,var_ANG_den,var_RANGE_den, varType, varNugget(:,3), [projectPath, inputPath, zoneFileName], nbrZones,...
%            nbrClassesBihist, 0, 'No File', 0, 'No file', 0, covTab),XX(1,1), YY(1,1), ZZ(1,1));

        % FT load LFM density as simulation  matrix
        [~,dgrid] = import_gslib([projectPath 'input\DeFlat.gslib']);
        DGrid = gslib2grid(dgrid,XX(1,1), YY(1,1), ZZ(1,1));
        
        twtsim=toc;
        fprintf(fileid, '%s %d \n', 'Simulations - elapsed time (sec)', twtsim); % WRITE NUMBER OF SIM
        clear twtsim

    else
% 
%         tic
%         fprintf(fileid, '%s %d \n', 'Simulating Density - Simulation number: ', noSim); % WRITE NUMBER OF SIM
%         DGrid = gslib2grid(parFileDSS(simType, 1, 1, 0, inputPath, den_hd, 1, [projectPath,outputPath, 'd_DSS_',num2str(noIt),'_',num2str(noSim)], 5,...
%            [projectPath, outputPath, 'best_D_',num2str(noIt-1),'.out'], repmat('No File', nbrZones,1), bounds_den,XX,YY,ZZ, ...
%            [projectPath,outputPath, 'best_cc_den',num2str(noIt-1),'.out'],0,'No File',krig_ANG_den,krig_RANGE_den,var_ANG_den,var_RANGE_den, varType, varNugget(:,1),...
%            [projectPath, inputPath, zoneFileName], nbrZones,...
%             nbrClassesBihist, 0, 'No File', 0, 'No file', 0, covTab),XX(1,1), YY(1,1), ZZ(1,1));
% 
%         fprintf(fileid, '%s %d \n', 'Simulating Vp - Simulation number: ', noSim); % WRITE NUMBER OF SIM       
%         VpGrid = gslib2grid(parFileDSS(simType, 1, 1, 1, inputPath, vp_hd, 1, [projectPath,outputPath, 'vp_coDSS_',num2str(noIt),'_',num2str(noSim)], 5,...
%            [projectPath,outputPath, 'best_vp_',num2str(noIt-1),'.out'], joint_d_vp, bounds_vp,XX,YY,ZZ, ...
%            [projectPath,outputPath, 'best_cc_vp',num2str(noIt-1),'.out'],0,[projectPath,outputPath, 'd_DSS_',num2str(noIt),'_',num2str(noSim),'_1.out']...
%            ,krig_ANG_vp,krig_RANGE_vp,var_ANG_vp,var_RANGE_vp, varType, varNugget(:,2), zoneFileName,nbrZones,...
%             nbrClassesBihist, 0, 'No File', 0, 'No file', 0, covTab),XX(1,1), YY(1,1), ZZ(1,1));
% 
%         fprintf(fileid, '%s %d \n', 'Simulating Vs - Simulation number: ', noSim); % WRITE NUMBER OF SIM
%         VsGrid = gslib2grid(parFileDSS(simType, 1, 1, 1, inputPath, vs_hd, 1, [projectPath,outputPath, 'vs_coDSS_',num2str(noIt),'_',num2str(noSim)], 5,...
%            [projectPath,outputPath, 'best_vs_',num2str(noIt-1),'.out'], joint_VsD , bounds_vs,XX,YY,ZZ, ...
%            [projectPath,outputPath, 'best_cc_vs',num2str(noIt-1),'.out'],0,[projectPath,outputPath, 'vp_coDSS_',num2str(noIt),'_',num2str(noSim),'_1.out']...
%            ,krig_ANG_vs,krig_RANGE_vs,var_ANG_vs,var_RANGE_vs, varType, varNugget(:,3), [projectPath, inputPath, zoneFileName],...
%            nbrZones, nbrClassesBihist, 0, 'No File', 0, 'No file', 0, covTab),XX(1,1), YY(1,1), ZZ(1,1));
% 
%         twtsim=toc;
%         fprintf(fileid, '%s %d \n', 'Simulations - elapsed time (sec)', twtsim); % WRITE NUMBER OF SIM
%         clear twtsim        
    end
    
    
end 
end